# 4-factor design with crossed random effects nitro <- read.csv('ecol 563/nitro.csv') library(lme4) mod3.lmer <- lmer(pN^2 ~ factor(lh)*factor(n)+factor(func)+factor(p)+(1|block)+(1|phy), data=nitro[nitro$tag!=444,]) # no p-values anova(mod3.lmer) # obtain credible intervals for effects out.mcmc <- mcmcsamp(mod3.lmer,n=10000) HPDinterval(out.mcmc) # effects model mod3.lm <- lm(pN^2 ~ factor(lh)*factor(n), data=nitro[nitro$tag!=444,]) summary(mod3.lm) # cell means model mod3.lm <- lm(pN^2 ~ factor(lh):factor(n)-1, data=nitro[nitro$tag!=444,]) # we get estimates and standard errors of means summary(mod3.lm) # variance-covariance matrix of means vcov(mod3.lm) # effect estimaters mod3.lm <- lm(pN^2 ~ factor(lh)+factor(n), data=nitro[nitro$tag!=444,]) summary(mod3.lm) # means for lh, effects for n mod3.lm <- lm(pN^2 ~ factor(lh)+factor(n)-1, data=nitro[nitro$tag!=444,]) summary(mod3.lm) # obtain means and their covariance matrix separately for different values of func and p fixef(mod3.lmer) # Method 1: create a matrix with 4 rows: one row for each lh and n combination cmat <- function(func,p) rbind(c(1,0,0,func=='mock',p==1,0), c(1,1,0,func=='mock',p==1,0), c(1,0,1,func=='mock',p==1,0), c(1,1,1,func=='mock',p==1,1)) cmat('inoc',0) # Method 2: write a function of all 4 factors that matches the R coding scheme cmat2 <- function(lh,n,func,p) c(1, lh=='NP', n==1, func=='mock', p==1, (lh=='NP')*(n==1)) # produce a matrix that contains all combinations of lh and n g <- expand.grid(lh=levels(nitro$lh), n=0:1) # evaluate function on this matrix apply(g,1, function(x) cmat2(x[1],x[2],'inoc',0)) # transpose result t(apply(g,1, function(x) cmat2(x[1],x[2],'inoc',0))) # write this as a function of func and p cmat3 <- function(func,p) t(apply(g,1,function(x) cmat2(x[1],x[2],func,p))) # we get the same result as above cmat3('inoc',0) #obtain estimates of the means ests <- cmat3('inoc',0) %*% fixef(mod3.lmer) #obtain the variance-covariance matrix of the means vmat <- cmat3('inoc',0) %*% vcov(mod3.lmer) %*% t(cmat3('inoc',0)) # modify the function for obtain difference-adjusted confidence intervals # change the t-functions to normal functions ci.func2 <- function(rowvals, lm.model, glm.vmat) { nor.func1a <- function(alpha, model, sig) 1-pnorm(-qnorm(1-alpha/2) * sum(sqrt(diag(sig))) / sqrt(c(1,-1) %*% sig %*%c (1,-1))) - pnorm(qnorm(1-alpha/2) * sum(sqrt(diag(sig))) / sqrt(c(1,-1) %*% sig %*% c(1,-1)), lower.tail=F) nor.func2 <- function(a,model,sigma) nor.func1a(a,model,sigma)-.95 n <- length(rowvals) xvec1b <- numeric(n*(n-1)/2) vmat <- glm.vmat[rowvals,rowvals] ind <- 1 for(i in 1:(n-1)) { for(j in (i+1):n){ sig <- vmat[c(i,j), c(i,j)] #solve for alpha xvec1b[ind] <- uniroot(function(x) nor.func2(x, lm.model, sig), c(.001,.999))$root ind <- ind+1 }} 1-xvec1b } # the variance-covariance matrix from lme4 has an unusual class class(vmat) # we need to change the class to type matrix class(as.matrix(vmat)) ci.func2(1:4, mod3.lmer, as.matrix(vmat)) # we get slightly diferent confidence levels depending on value of func and p vmat2 <- cmat3('mock',0) %*% vcov(mod3.lmer) %*% t(cmat3('mock',0)) ci.func2(1:4, mod3.lmer, as.matrix(vmat2)) vmat2 <- cmat3('mock',1) %*% vcov(mod3.lmer) %*% t(cmat3('mock',1)) ci.func2(1:4,mod3.lmer,as.matrix(vmat2)) # function to obtain means, std errs, 95% CIs, and difference-adjusted CIs make.data <- function(func,p,model) { vmat <- cmat3(func,p) %*% vcov(model) %*% t(cmat3(func,p)) vals1 <- ci.func2(1:4, model, as.matrix(vmat)) part1 <- data.frame(lh=c('EA','NP','EA','NP'), n=c(0,0,1,1), est=as.vector(cmat3(func,p)%*%fixef(model)), se=sqrt(diag(vmat))) part1$low95 <- part1$est+qnorm(.025)*part1$se part1$up95 <- part1$est+qnorm(.975)*part1$se part1$low85 <- part1$est+qnorm((1-min(vals1))/2)*part1$se part1$up85 <- part1$est+qnorm(1-(1-min(vals1))/2)*part1$se part1$func <- func part1$p <- p part1 } # use function for different values of func and p part0a <- make.data('inoc',0,mod3.lmer) part0b <- make.data('inoc',1,mod3.lmer) part0c <- make.data('mock',0,mod3.lmer) part0d <- make.data('mock',1,mod3.lmer) # assemble results in single matrix fac.vals <- rbind(part0a,part0b,part0c,part0d) fac.vals # there is no function for getting confidence intervals for lme4 objects # we can obtain confidence intervals for an lme object with the intervals function library(nlme) mod3.lme <- lme(pN^2 ~ factor(lh)*factor(n)+factor(func)+factor(p), random=~1|block, data=nitro[nitro$tag!=444,]) # confint does not work on lme objects confint(mod3.lme) methods(confint) # intervals works intervals(mod3.lme) # if we had a simpler model we could use intervals to obtain confidence intervals for mean # effects version of 2-factor interaction model mod3.lme <- lme(pN^2 ~ factor(lh)*factor(n),random=~1|block, data=nitro[nitro$tag!=444,]) # means version of 2-factor interaction model mod3a.lme <- lme(pN^2 ~ factor(lh):factor(n)-1,random=~1|block, data=nitro[nitro$tag!=444,]) # obtain means fixef(mod3a.lme) # obtain confidence intervals intervals(mod3a.lme) # obtain variance covariance matrix of means vcov(mod3a.lme) # group panel function my.panel <- function(x, y, subscripts, col, pch, group.number, ...) { #subset variables for the current panel low95 <- fac.vals$low95[subscripts] up95 <- fac.vals$up95[subscripts] low85 <- fac.vals$low85[subscripts] up85 <- fac.vals$up85[subscripts] myjitter <- c(-.05,.05) col2 <- c('tomato','grey60') #95% confidence interval panel.arrows(x+myjitter[group.number],low95,x+myjitter[group.number],up95,angle=90, code=3, length=.05, col=col) #50% confidence interval panel.segments(x+myjitter[group.number],low85,x+myjitter[group.number],up85, col=col2[group.number], lineend=1, lwd=5) panel.lines(x+myjitter[group.number],y, col=col) panel.xyplot(x+myjitter[group.number],y, col='white', pch=16) panel.xyplot(x+myjitter[group.number],y, col=col, pch=pch, ...) panel.abline(h=seq(10,26,2), col='lightgrey', lty=3) } # prepanel function to set limits for panels prepanel.ci2 <- function(x,y,ly,uy,subscripts,...){ list(ylim=range(y, uy[subscripts], ly[subscripts], finite=T))} library(lattice) # two-factor interaction plot using lattice and group panel function xyplot(est~lh|func+factor(p, level=0:1, labels=paste('P = ', 0:1, sep='')), xlab ='lh', ylab='Estimated mean', groups=factor(n), data=fac.vals, col=c(2,1), pch=c(1,16), prepanel=prepanel.ci2, ly=fac.vals$low95, uy=fac.vals$up95, layout=c(2,2), panel.groups=my.panel, panel="panel.superpose") # add a key to the graph xyplot(est~lh|func+factor(p, level=0:1, labels=paste('P = ', 0:1, sep='')), xlab ='lh', ylab='Estimated mean', groups=factor(n), data=fac.vals, col=c(2,1), pch=c(1,16), prepanel=prepanel.ci2, ly=fac.vals$low95, uy=fac.vals$up95, layout=c(2,2), panel.groups=my.panel, panel="panel.superpose", key=list(x=.85,y=.80, corner=c(0,0), points=list(pch=c(16,1), col=1:2, cex=.9), text=list(c('n = 1','n = 0'), cex=.8))) # assign the graph to an object mygraph <- xyplot(est~lh|func+factor(p, level=0:1, labels=paste('P = ', 0:1, sep='')), xlab ='lh', ylab='Estimated mean', groups=factor(n), data=fac.vals, col=c(2,1), pch=c(1,16), prepanel=prepanel.ci2, ly=fac.vals$low95, uy=fac.vals$up95, layout=c(2,2), panel.groups=my.panel, panel="panel.superpose", key=list(x=.85,y=.88, corner=c(0,0), points=list(pch=c(16,1), col=1:2, cex=.9), text=list(c('n = 1','n = 0'), cex=.8))) # load the latticeExtra package and redo graph library(latticeExtra) # put strips on left and on top useOuterStrips(mygraph)